###########################################################################
###########################################################################
####Non-state Atrocities in Capital Cities - ZINB Urban Sample Analysis####
###########################################################################
###########################################################################
library(MASS) 
library(pscl)
library(foreign)
library(stargazer)
library(mvtnorm)
library(plyr)

##################################
###Create Urban Cells Only Data###
##################################
# Set working library
setwd("~/Data/Global Analysis/")
# Read in main data
main.data <- read.dta("full.grid.upd.dta")
#subset out only cells with some urbanization levels
u.main.data <- main.data[ which(main.data$lagurban > 0),]
#Identify threshold for large urban area (95 percentile)
quantile(u.main.data$lagurban, 0.95)
#Create an indicator for large urban areas
u.main.data$lag.large.urb <- ifelse(u.main.data$lagurban>=quantile(u.main.data$lagurban, 0.95),1,0)

##################################
###Global Urban Sample Analysis###
##################################
###Model XIII-A (large urban areas only in count stage)
at.zinb.ai.xiii.a <- zeroinfl(incidentnonstatefull ~ lag_Capital + lagcivconflagtemp + lagincidentsumfull + lag.large.urb + loglagross_oil_prod
                        + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + year96 + year97+
                          year98 + year99 + year00 + year01 + year02 + year03 + year04 + year05 + year06 + year07 + year08
                        |lagcivconflagtemp + lagincidentsumfull + loglagross_oil_prod
                        + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea, 
                        dist = "negbin", data = u.main.data)
summary(at.zinb.ai.xiii.a)
AIC(at.zinb.ai.xiii.a)

###Model XIII-B (large urban areas in both stages)
at.zinb.ai.xiii.b <- zeroinfl(incidentnonstatefull ~ lag_Capital + lagcivconflagtemp + lagincidentsumfull + lag.large.urb + loglagross_oil_prod
                        + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + year96 + year97+
                          year98 + year99 + year00 + year01 + year02 + year03 + year04 + year05 + year06 + year07 + year08
                        |lagcivconflagtemp + lagincidentsumfull + lag.large.urb + loglagross_oil_prod
                        + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea, 
                        dist = "negbin", data = u.main.data)
summary(at.zinb.ai.xiii.b)
AIC(at.zinb.ai.xiii.b)

###Model XIII-C (With urban in inflation stage)
at.zinb.ai.xiii.c <- zeroinfl(incidentnonstatefull ~ lag_Capital + lagcivconflagtemp + lagincidentsumfull + lag.large.urb + loglagross_oil_prod
                        + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + year96 + year97+
                          year98 + year99 + year00 + year01 + year02 + year03 + year04 + year05 + year06 + year07 + year08
                        |lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                        + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea, 
                        dist = "negbin", data = u.main.data)
summary(at.zinb.ai.xiii.c)
AIC(at.zinb.ai.xiii.c)

##################################
###Africa Urban Sample Analysis###
##################################
###Subset out African country cells only
u.main.data.Af<-subset(u.main.data,(u.main.data$ccode>400 & u.main.data$ccode<630) | u.main.data$ccode==651)

###Model XIV-A (large urban areas only in count stage, corresponding to Model 5A)
at.zinb.ai.xiv.a <- zeroinfl(incidentnonstatefull ~ lag_Capital + lagcivconflagtemp + lagincidentsumfull + lag.large.urb + loglagross_oil_prod
                        + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + year96 + year97+
                          year98 + year99 + year00 + year01 + year02 + year03 + year04 + year05 + year06 + year07 + year08
                        |lagcivconflagtemp + lagincidentsumfull + loglagross_oil_prod
                        + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea, 
                        dist = "negbin", data = u.main.data.Af)
summary(at.zinb.ai.xiv.a)
AIC(at.zinb.ai.xiv.a)

###Model XIV-B (large urban areas in both stages, corresponding to Model 5B)
at.zinb.ai.xiv.b <- zeroinfl(incidentnonstatefull ~ lag_Capital + lagcivconflagtemp + lagincidentsumfull + lag.large.urb + loglagross_oil_prod
                        + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + year96 + year97+
                          year98 + year99 + year00 + year01 + year02 + year03 + year04 + year05 + year06 + year07 + year08
                        |lagcivconflagtemp + lagincidentsumfull + lag.large.urb + loglagross_oil_prod
                        + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea, 
                        dist = "negbin", data = u.main.data.Af)
summary(at.zinb.ai.xiv.b)
AIC(at.zinb.ai.xiv.b)

###Model XIV-C (With urban in inflation stage, corresponding to Model 5C)
at.zinb.ai.xiv.c <- zeroinfl(incidentnonstatefull ~ lag_Capital + lagcivconflagtemp + lagincidentsumfull + lag.large.urb + loglagross_oil_prod
                        + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + year96 + year97+
                          year98 + year99 + year00 + year01 + year02 + year03 + year04 + year05 + year06 + year07 + year08
                        |lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                        + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea, 
                        dist = "negbin", data = u.main.data.Af)
summary(at.zinb.ai.xiv.c)
AIC(at.zinb.ai.xiv.c)

#########Export all models to LaTex
stargazer(at.zinb.ai.xiii.a, at.zinb.ai.xiii.b, at.zinb.ai.xiii.c, at.zinb.ai.xiv.a, at.zinb.ai.xiv.b, at.zinb.ai.xiv.c)
stargazer(at.zinb.ai.xiii.a, at.zinb.ai.xiii.b, at.zinb.ai.xiii.c, at.zinb.ai.xiv.a, at.zinb.ai.xiv.b, at.zinb.ai.xiv.c, zero.component=TRUE)
AIC(at.zinb.ai.xiii.a, at.zinb.ai.xiii.b, at.zinb.ai.xiii.c, at.zinb.ai.xiv.a, at.zinb.ai.xiv.b, at.zinb.ai.xiv.c)

##################################################
####Vuong Tests -- Negative Binomial Comparison###
##################################################
##Models XIII-A -- XIII-C
nb.xiii <- glm.nb(incidentnonstatefull ~ lag_Capital + lagcivconflagtemp + lagincidentsumfull + lag.large.urb + loglagross_oil_prod
              + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + year96 + year97+
                year98 + year99 + year00 + year01 + year02 + year03 + year04 + year05 + year06 + year07 + year08, 
              data = u.main.data)
#Compare NB and ZINB
vuong(nb.xiii, at.zinb.ai.xiii.a)
vuong(nb.xiii, at.zinb.ai.xiii.b)
vuong(nb.xiii, at.zinb.ai.xiii.c)

##Models XIV-A -- XIV-C
nb.xiv <- glm.nb(incidentnonstatefull ~ lag_Capital + lagcivconflagtemp + lagincidentsumfull + lag.large.urb + loglagross_oil_prod
              + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + year96 + year97+
                year98 + year99 + year00 + year01 + year02 + year03 + year04 + year05 + year06 + year07 + year08, 
              data = u.main.data.Af)
#Compare NB and ZINB
vuong(nb.xiv, at.zinb.ai.xiv.a)
vuong(nb.xiv, at.zinb.ai.xiv.b)
vuong(nb.xiv, at.zinb.ai.xiv.c)
